## Import data from .RData
load("data/cont_snotel.RData")
# Import seasonal snow zone polygons
SSZ <- st_read('data/MODIS_Snow_Zones/SSZ_0cc.shp')
## Reading layer `SSZ_0cc' from data source
## `C:\Users\wreis\OneDrive - Colostate\MS Research\R_Git\SNOTEL\data\MODIS_Snow_Zones\SSZ_0cc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2357223 ymin: 1163534 xmax: -437223.4 ymax: 3151034
## Projected CRS: NAD_1983_Albers
PSZ <- st_read('data/MODIS_Snow_Zones/PSZ_0cc.shp')
## Reading layer `PSZ_0cc' from data source
## `C:\Users\wreis\OneDrive - Colostate\MS Research\R_Git\SNOTEL\data\MODIS_Snow_Zones\PSZ_0cc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2357223 ymin: 1459534 xmax: -437223.4 ymax: 3149702
## Projected CRS: NAD_1983_Albers
# Group sites by site_id
cont_snotel_site <- cont_snow_data %>%
group_by(site_id, site_name, state, latitude, longitude, elev) %>%
summarise()
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# transform files to CRS:4326
snotel_sites <- st_as_sf(cont_snotel_site, coords = c('longitude', 'latitude'), crs = 4326) #fix before next run
PSZ_4326 <- st_transform(PSZ, crs = 4326)
# clip SNOTEL Sites to PSZ
sf::sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
PSZ_snotel_sites = st_intersection(PSZ_4326, snotel_sites)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
mapview(PSZ_snotel_sites) +
mapview(PSZ_4326)
# Filter all SNOTEL to PSZ SNOTEL Sites Only
PSZ_cont_snotel <- cont_snow_data[cont_snow_data$site_id %in% PSZ_snotel_sites$site_id,]
# PSZ SNOTEL Site, Calculate spring days with new SWE and new SWE depth
PSZ_cont_snotel_spring <- PSZ_cont_snotel %>%
mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
filter(state != "SD", year >= 1979, ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244))
# SNOTEL Analysis, calculate the number and amount of positive SWE between
cont_snotel_spring <- cont_snow_data %>%
mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
filter(state != "SD", year >= 1979, ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244))
# PSZ determine average number of days with increased SWE per spring at each site
PSZ_cont_snotel_site <- PSZ_cont_snotel_spring %>%
group_by(site_id, site_name, state, latitude, longitude, elev) %>%
summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), years = length(unique(year)), avg_spring_days = newsnow_days/years,
newsnow_depth = sum(newsnow, na.rm = TRUE), avg_spring_snow = newsnow_depth/years)
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Determine average number of days with increased SWE per spring at each site
cont_snotel_site <- cont_snotel_spring %>%
group_by(site_id, site_name, state, latitude, longitude, elev) %>%
summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), years = length(unique(year)), avg_spring_days = newsnow_days/years,
newsnow_depth = sum(newsnow, na.rm = TRUE), avg_spring_snow = newsnow_depth/years)
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# View data
mapview(cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_days",
layer.name = "AVG Sping Days (All Sites)" ,crs = 4326, grid = FALSE) +
mapview(PSZ_cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_days",
layer.name = "AVG Sping Days (PSZ Sites)" , crs = 4326, grid = FALSE)
# Days of new snow at SNOTEL Sites in Western US
#Summarize the number of new SWE days and SWE depth at all sites per year
cont_snotel_site_yr <- cont_snotel_spring %>%
group_by(site_id, site_name, state, year, latitude, longitude, elev) %>%
summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), newsnow_depth = sum(newsnow, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Create a box plot for all sites
spring_snow_days <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_days, color = state)) +
geom_boxplot()+
ggtitle("Days of increased SWE at all sites (1978-2022)")
ggplotly(spring_snow_days)
#Summarize the number of new SWE days and SWE depth at PSZ sites per year
PSZ_cont_snotel_site_yr <- PSZ_cont_snotel_spring %>%
group_by(site_id, site_name, state, year, latitude, longitude, elev) %>%
summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), newsnow_depth = sum(newsnow, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Create a boxplot for all sites
PSZ_spring_snow_days <- ggplot(PSZ_cont_snotel_site_yr, aes(x = state, y= newsnow_days, color = state)) +
geom_boxplot()+
ggtitle("Days of increased SWE at PSZ sites (1978-2022)")
ggplotly(PSZ_spring_snow_days)
summary(cont_snotel_site_yr$elev)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 128 1817 2393 2290 2772 3544
hist(cont_snotel_site_yr$elev)

cont_snotel_site_yr <- cont_snotel_site_yr %>%
mutate(elev_band = ifelse(elev >= 2000, "High (>=2500)",
ifelse(elev < 2000 & elev >= 1200, "Mid (1500-2500)", "Low (<= 1500)")))
# Elevation vs new snow count
elev_newsnow <- plot_ly(cont_snotel_site, x = ~avg_spring_days, y = ~elev, color = ~state, type = 'scatter', mode = 'markers',
hoverinfo = 'text',
text = ~paste('</br> avg_new_days: ', avg_spring_days,
'</br> elev: ', elev,
'</br> site_name: ', site_name,
'</br> state: ', state)) %>%
layout(title = "Average Increased SWE Days by Elevation")
elev_newsnow
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
elev_newsnow_depth <- plot_ly(cont_snotel_site, x = ~avg_spring_snow, y = ~elev, color = ~state, type = 'scatter', mode = 'markers',
hoverinfo = 'text',
text = ~paste('</br> avg_new_SWE (mm): ', avg_spring_snow,
'</br> elev: ', elev,
'</br> site_name: ', site_name,
'</br> state: ', state)) %>%
layout(title = "Average Increased SWE Depth (mm) by Elevation")
elev_newsnow_depth
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
elev_boxplot <- ggplot(cont_snotel_site_yr, aes(x = elev_band, y= newsnow_days, color = elev_band)) +
geom_boxplot()+
ggtitle("New Snow Events by Elevation")
ggplotly(elev_boxplot)
# Spring Snow depth
spring_snow_depth <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_depth, color = state)) +
geom_boxplot()+
ggtitle("Avg Depth of increased SWE at all sites (1978-2022)")
ggplotly(spring_snow_depth)
# Spring Snow depth
PSZ_spring_snow_depth <- ggplot(PSZ_cont_snotel_site_yr, aes(x = state, y= newsnow_depth, color = state)) +
geom_boxplot()+
ggtitle("Depth of increased SWE at PSZ sites (1978-2022)")
ggplotly(PSZ_spring_snow_depth)
# spring_snow_days_yr <- ggplot(cont_snotel_site_yr, aes(x = year, y = newsnow_days, group = state, color = state)) +
# geom_line()
#
# ggplotly(spring_snow_days_yr)